import numpy as np
from spatio_flux import SpatioFluxVivarium, build_path, render_path
Make a Vivarium¶
view available types and processes
# TODO -- make a spatio-flux specific Vivarium with the the core preloaded
vi = SpatioFluxVivarium()
# view the available types
vi.get_types()
['', 'string', 'length*mass/current^2*time^2', 'mass/time^3', 'list', 'process', 'length^2*mass/time^3', 'reaction', 'step', 'bounds', 'mass/current*time^2', 'length^2*mass/temperature*time^2', 'length^2*mass/time^2', 'printing_unit', 'length/time', 'time', 'length*mass/time^2', 'current*time/substance', 'length^4*mass/current*time^3', 'emitter_mode', 'number', 'time^2/length', 'kinetics', '/substance', 'length^2', 'current^2*time^4/length^2*mass', 'path', 'current*time', 'current*length^2*time', 'current*length^2', 'temperature', 'luminosity/length^2', 'interval', 'current^2*time^4/length^3*mass', 'substrate_role', 'length^0_5*mass^0_5/time', 'length^1_5*mass^0_5/time', '/temperature*time', 'tree', 'wires', 'length^2*mass/time', 'quote', 'current*length*time', 'length/mass', 'current*time/mass', 'substance', 'map', '/printing_unit', 'tuple', 'mass/length*time', 'length^2*mass/current^2*time^2', '/time', 'mass', 'edge', 'boundary_side', 'current', 'float', 'length^3', 'length^2*mass/current*time^3', 'mass^0_5/length^1_5', 'maybe', 'enum', 'positive_float', 'length^3*mass/current^2*time^4', 'length^2*mass/substance*temperature*time^2', 'length*temperature', 'length^1_5*mass^0_5/time^2', 'luminosity', 'length^2/time^2', 'mass/temperature^4*time^3', 'mass/length', 'protocol', 'time/length', 'union', 'schema', 'current^2*time^3/length^2*mass', 'boolean', 'length*time/mass', 'integer', 'composite', 'length', 'length^2*mass/current^2*time^3', 'length^4*mass/time^3', 'array', 'length^2/time', 'substance/length^3', 'printing_unit/length', '/length', 'length^0_5*mass^0_5', 'any', 'length^3/mass*time^2', 'current*time^2/length^2*mass', 'mass^0_5/length^0_5*time', 'length^3/time', 'substance/time', 'length^2*mass/current*time^2', 'mass/length*time^2', 'length*mass/current*time^3', 'mass/length^3', 'mass/time^2', 'particle', 'length/time^2']
# view the available processes
vi.get_processes()
['composite', 'ram-emitter', 'DiffusionAdvection', 'console-emitter', 'DynamicFBA', 'MinimalParticle', 'Particles', 'json-emitter']
dFBA¶
Dynamic Flux Balance Analysis (dFBA) extends traditional Flux Balance Analysis (FBA) to model the dynamic behavior of metabolic networks over time, allowing for the simulation of growth and substrate consumption in a changing environment.
# inspect the config schema for the 'increase' process
vi.process_config('DynamicFBA')
{'model_file': 'string',
'kinetic_params': 'map[tuple[float,float]]',
'substrate_update_reactions': 'map[string]',
'biomass_identifier': 'string',
'bounds': 'map[bounds]'}
vi.core.default(vi.process_config('DynamicFBA')) # TODO use default to
{'model_file': '',
'kinetic_params': {},
'substrate_update_reactions': {},
'biomass_identifier': '',
'bounds': {}}
# dfba_config = vi.process_config('DynamicFBA', dataclass=True) # TODO -- make this work
# TODO - enable get_dataclass to work with the new process
# dfba_config = v.get_dataclass('DynamicFBA')
dfba_config = {
"model_file": "textbook",
"kinetic_params": {
"glucose": (0.5, 1),
"acetate": (0.5, 2)},
"substrate_update_reactions": {
"glucose": "EX_glc__D_e",
"acetate": "EX_ac_e"},
"biomass_identifier": "biomass",
"bounds": {
"EX_o2_e": {"lower": -2, "upper": None},
"ATPM": {"lower": 1, "upper": 1}
}
}
# make the vivarium
v = SpatioFluxVivarium()
# add a dynamic FBA process called 'dFBA'
v.add_process(name="dFBA",
process_id="DynamicFBA",
config=dfba_config,
)
v.diagram(dpi='70')
mol_ids = ["glucose", "acetate", "biomass"]
path=["fields"]
for mol_id in mol_ids:
v.add_object(
name=mol_id,
path=path,
value=np.random.rand()
)
v.connect_process(
name="dFBA",
inputs={
"substrates": {
mol_id: build_path(path, mol_id)
for mol_id in mol_ids}
},
outputs={
"substrates": {
mol_id: build_path(path, mol_id)
for mol_id in mol_ids}
}
)
# add an emitter to save results
v.add_emitter()
# TODO -- show_values does not work
v.diagram(dpi='70', show_values=True)
v.set_value(path = ['fields', 'glucose'], value=10)
v.set_value(path = ['fields', 'biomass'], value=0.1)
field = v.get_value(['fields'])
print(field)
{'glucose': 10, 'acetate': 0.6199272885535164, 'biomass': 0.1}
v.save(filename='dFBA_t0')
Saved file: out/dFBA_t0.json
# run the simulation for 10 time units
v.run(interval=60)
v.plot_timeseries(
subplot_size=(8, 3),
combined_vars=[
[
'/fields/glucose',
'/fields/acetate',
'/fields/biomass'
]
]
)
Spatial dFBA¶
mol_ids = ["glucose", "acetate", "biomass"]
path=["fields"]
rows = 2
columns = 1
# make the vivarium
v2 = SpatioFluxVivarium()
for mol_id in mol_ids:
v2.add_object(
name=mol_id,
path=path,
value=np.random.rand(rows, columns)
)
# add a dynamic FBA process at every location
for i in range(rows):
for j in range(columns):
dfba_name = f"dFBA[{i},{j}]"
# add a process for this location
v2.add_process(
name=dfba_name,
process_id="DynamicFBA",
config=dfba_config,
)
v2.connect_process(
name=dfba_name,
inputs={"substrates": {
mol_id: build_path(path, mol_id, i, j)
for mol_id in mol_ids}},
outputs={"substrates": {
mol_id: build_path(path, mol_id, i, j)
for mol_id in mol_ids}}
)
# add an emitter to save results
v2.add_emitter()
v2.diagram(dpi='70')
# change some initial values
v2.merge_value(path = ['fields', 'glucose', 0, 0], value=10.0)
v2.merge_value(path = ['fields', 'biomass', 0, 0], value=0.1)
field = v2.get_value(['fields'])
print(field)
{'glucose': array([[10. ],
[ 0.13101474]]), 'acetate': array([[0.47702186],
[0.37300945]]), 'biomass': array([[0.1 ],
[0.56016474]])}
v2.run(60)
v2.get_timeseries(as_dataframe=True)
| /global_time | /fields/glucose | /fields/acetate | /fields/biomass | |
|---|---|---|---|---|
| 0 | 0.0 | [[10.0], [0.1310147415999353]] | [[0.47702185787766227], [0.37300945027621435]] | [[0.1], [0.5601647442497885]] |
| 1 | 1.0 | [[9.904761904761905], [0.014710265043333406]] | [[0.4891299100737515], [0.0]] | [[0.1079432568708625], [0.578804696512899]] |
| 2 | 2.0 | [[9.80200585245463], [0.014710265043333406]] | [[0.5021057975498924], [0.0]] | [[0.11651530697082768], [0.578804696512899]] |
| 3 | 3.0 | [[9.691145527078628], [0.014710265043333406]] | [[0.5160006282993799], [0.0]] | [[0.12576552148631978], [0.578804696512899]] |
| 4 | 4.0 | [[9.571550338512791], [0.014710265043333406]] | [[0.530866005462576], [0.0]] | [[0.13574706721093277], [0.578804696512899]] |
| ... | ... | ... | ... | ... |
| 56 | 56.0 | [[0.0], [0.014710265043333406]] | [[0.0], [0.0]] | [[0.9748198839578434], [0.578804696512899]] |
| 57 | 57.0 | [[0.0], [0.014710265043333406]] | [[0.0], [0.0]] | [[0.9748198839578434], [0.578804696512899]] |
| 58 | 58.0 | [[0.0], [0.014710265043333406]] | [[0.0], [0.0]] | [[0.9748198839578434], [0.578804696512899]] |
| 59 | 59.0 | [[0.0], [0.014710265043333406]] | [[0.0], [0.0]] | [[0.9748198839578434], [0.578804696512899]] |
| 60 | 60.0 | [[0.0], [0.014710265043333406]] | [[0.0], [0.0]] | [[0.9748198839578434], [0.578804696512899]] |
61 rows × 4 columns
# get a list of all the paths so they can be plotted together in a single graph
all_paths = []
for i in range(rows):
for j in range(columns):
# get the paths for this location
location_path = []
for mol_id in mol_ids:
this_path = build_path(path, mol_id, i, j)
rendered_path = render_path(this_path)
location_path.append(rendered_path)
all_paths.append(location_path)
# print(all_paths)
v2.plot_timeseries(
subplot_size=(8, 3),
combined_vars=all_paths
)
v2.show_video()
Diffusion/Advection¶
This approach models the physical processes of diffusion and advection in two dimensions, providing a way to simulate how substances spread and are transported across a spatial domain, essential for understanding patterns of concentration over time and space.
vi.core.default(vi.process_config('DiffusionAdvection'))
{'n_bins': (0, 0),
'bounds': (0.0, 0.0),
'default_diffusion_rate': 0.1,
'default_diffusion_dt': 0.1,
'diffusion_coeffs': {},
'advection_coeffs': {}}
bounds = (10.0, 10.0)
n_bins = (10, 10)
mol_ids = [
'glucose',
'acetate',
'biomass'
]
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {
'biomass': (0, -0.1)
}
path = ['fields']
v3 = SpatioFluxVivarium()
for mol_id in mol_ids:
v3.add_object(
name=mol_id,
path=path,
value=np.random.rand(n_bins[0], n_bins[1])
)
v3.add_process(name='diffusion_advection',
process_id='DiffusionAdvection',
config={
'n_bins': n_bins,
'bounds': bounds,
'default_diffusion_rate': diffusion_rate,
'default_diffusion_dt': diffusion_dt,
# 'diffusion_coeffs': diffusion_coeffs_all,
'advection_coeffs': advection_coeffs,
},
inputs={'fields': ['fields']},
outputs={'fields': ['fields']}
)
# add an emitter to save results
v3.add_emitter()
v3.diagram(dpi='70')
v3.run(60)
v3.plot_snapshots(n_snapshots=5)
v3.show_video()
COMETS (Computation Of Microbial Ecosystems in Time and Space)¶
COMETS combines dynamic FBA with spatially resolved physical processes (like diffusion and advection) to simulate the growth, metabolism, and interaction of microbial communities within a structured two-dimensional environment, capturing both biological and physical complexities.
bounds = (20.0, 10.0) # Bounds of the environment
n_bins = (12, 6)
mol_ids = ['glucose', 'acetate', 'biomass']
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {
'biomass': (0, 0.1)
}
path = ['fields']
# Initialize the acetate concentration field with zeros
acetate_field = np.zeros((n_bins[0], n_bins[1]))
# Create a linear glucose concentration gradient from 1 (top) to 0 (bottom) vertically
max_glc = 10
glc_field = np.random.rand(n_bins[0], n_bins[1]) * max_glc
# Initialize the biomass concentration field, setting a small biomass concentration in a specific row
biomass_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field[0:int(1*n_bins[0]/5), int(2*n_bins[1]/5):int(3*n_bins[1]/5)] = 0.1 # place some biomass
# make the vivarium
v4 = SpatioFluxVivarium()
# create the molecular fields
# for mol_id in mol_ids:
v4.add_object(
name='glucose',
path=path,
value=glc_field
)
v4.add_object(
name='biomass',
path=path,
value=biomass_field
)
v4.add_object(
name='acetate',
path=path,
value=acetate_field
)
# add a diffusion/advection process
v4.add_process(name='diffusion_advection',
process_id='DiffusionAdvection',
config={
'n_bins': n_bins,
'bounds': bounds,
'default_diffusion_rate': diffusion_rate,
'default_diffusion_dt': diffusion_dt,
'advection_coeffs': advection_coeffs},
inputs={'fields': path},
outputs={'fields': path}
)
# add a dynamic FBA process at every location
for i in range(n_bins[0]):
for j in range(n_bins[1]):
dfba_name = f"dFBA[{i},{j}]"
# add a process for this location
v4.add_process(
name=dfba_name,
process_id="DynamicFBA",
config=dfba_config,
)
v4.connect_process(
name=dfba_name,
inputs={"substrates": {
mol_id: build_path(path, mol_id, i, j)
for mol_id in mol_ids}},
outputs={"substrates": {
mol_id: build_path(path, mol_id, i, j)
for mol_id in mol_ids}}
)
# add an emitter to save results
v4.add_emitter()
v4.diagram(dpi='70',
remove_nodes=[f"/dFBA[{i},{j}]" for i in range(n_bins[0]-1) for j in range(n_bins[1])]
)
v4.run(60)
v4.plot_snapshots(n_snapshots=5)
v4.show_video()
v4.plot_timeseries(
subplot_size=(8, 3),
query=[
'/fields/glucose/0/0',
'/fields/acetate/0/0',
'/fields/biomass/0/0',
],
combined_vars=[[
'/fields/glucose/0/0',
'/fields/acetate/0/0',
'/fields/biomass/0/0',
]]
)
Particles¶
# TODO -- remove this from final
import numpy as np
from spatio_flux import SpatioFluxVivarium
bounds = (10.0, 20.0) # Bounds of the environment
n_bins = (20, 40) # Number of bins in the x and y directions
v5 = SpatioFluxVivarium()
v5.add_process(
name='particle_movement',
process_id='Particles',
config={
'n_bins': n_bins,
'bounds': bounds,
'diffusion_rate': 0.1,
'advection_rate': (0, -0.1),
'add_probability': 0.4,
'boundary_to_add': ['top']
},
)
v5.connect_process(
name='particle_movement',
inputs={
'fields': ['fields'],
'particles': ['particles']
},
outputs={
'fields': ['fields'],
'particles': ['particles']
}
)
v5.initialize_process(
path='particle_movement',
config={'n_particles': 2}
)
v5.add_emitter()
v5.diagram(dpi='70')
v5.run(100)
v5_results = v5.get_results()
from spatio_flux.viz.plot import plot_species_distributions_with_particles_to_gif
# TODO -- integrate this method with vivarium
plot_species_distributions_with_particles_to_gif(
v5_results,
skip_frames=3,
bounds=bounds
)
Saving GIF to species_distribution_with_particles.gif
v5.diagram(dpi='70')
v5.save(filename='v5_post_run', outdir='out')
Saved file: out/v5_post_run.json
Minimal particle process¶
from spatio_flux.processes.particles import get_minimal_particle_composition
particle_schema = get_minimal_particle_composition(v5.core)
document = v5.make_document()
# document['state']['particles'] = {}
# # document['composition']['particles']
# document['composition'] = particle_schema
document['state']['global_time'] = 0.0
v6 = SpatioFluxVivarium(document=document)
v6.diagram(dpi='70')
v6.run(10)
v6_results = v6.get_results()
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) Cell In[87], line 1 ----> 1 v6.run(10) 2 v6_results = v6.get_results() File ~/code/vivarium-interface/vivarium/vivarium.py:564, in Vivarium.run(self, interval) 561 if not self.emitter_paths: 562 self.add_emitter() --> 564 self.composite.run(interval) File ~/code/process-bigraph/process_bigraph/composite.py:783, in Composite.run(self, interval, force_complete) 781 for path in self.process_paths: 782 process = get_path(self.state, path) --> 783 full_step = self.run_process( 784 path, 785 process, 786 end_time, 787 full_step, 788 force_complete) 790 # apply updates based on process times in self.front 791 if full_step == math.inf: 792 # no processes ran, jump to next process File ~/code/process-bigraph/process_bigraph/composite.py:680, in Composite.run_process(self, path, process, end_time, full_step, force_complete) 677 full_step = interval 679 if future <= end_time: --> 680 update = self.process_update( 681 path, 682 process, 683 state, 684 process_interval 685 ) 687 # update front, to be applied at its projected time 688 self.front[path]['time'] = future File ~/code/process-bigraph/process_bigraph/composite.py:622, in Composite.process_update(self, path, process, states, interval, ports_key) 600 """Start generating a process's update. 601 602 This function is similar to :py:meth:`_invoke_process` except in (...) 617 ``store``. 618 """ 620 states = strip_schema_keys(states) --> 622 update = process['instance'].invoke(states, interval) 624 def defer_project(update, args): 625 schema, state, path = args File ~/code/process-bigraph/process_bigraph/composite.py:303, in Process.invoke(self, state, interval) 302 def invoke(self, state, interval): --> 303 update = self.update(state, interval) 304 sync = SyncUpdate(update) 305 return sync File ~/code/spatio-flux/spatio_flux/processes/particles.py:183, in Particles.update(self, state, interval) 180 # Apply diffusion and advection 181 dx, dy = np.random.normal(0, self.config['diffusion_rate'], 2) + self.config['advection_rate'] --> 183 new_x_position = particle['position'][0] + dx 184 new_y_position = particle['position'][1] + dy 186 # Check and remove particles if they hit specified boundaries TypeError: can only concatenate str (not "numpy.float64") to str
from spatio_flux.processes.particles import get_minimal_particle_composition
bounds = (10.0, 20.0) # Bounds of the environment
n_bins = (4, 8) # Number of bins in the x and y directions
# same as before
v6 = SpatioFluxVivarium()
v6.add_process(
name='particle_movement',
process_id='Particles',
config={
'n_bins': n_bins,
'bounds': bounds,
'diffusion_rate': 0.1,
'advection_rate': (0, -0.1),
'add_probability': 0.4,
'boundary_to_add': ['top']})
v6.connect_process(
name='particle_movement',
inputs={
'fields': ['fields'],
'particles': ['particles']},
outputs={
'fields': ['fields'],
'particles': ['particles']})
# Here we go beyond the previous particle simulation
# add a minimal particle process into the schema
particle_schema = get_minimal_particle_composition(v6.core)
v6.set_schema('particles', particle_schema['particles'])
v6.initialize_process(
path='particle_movement',
config={'n_particles': 1})
v6.diagram(dpi='70')
# v6
# particle_schema = get_minimal_particle_composition(v6.core)
# particle_schema
v6.run(10)
v6_results = v6.get_results()
from spatio_flux.viz.plot import plot_species_distributions_with_particles_to_gif
# TODO -- integrate this method with vivarium
plot_species_distributions_with_particles_to_gif(
v6_results,
skip_frames=2,
bounds=bounds
)